Constitutive Modeling with Single and Dual Internal Variables

Phenomenological constitutive models with internal variables have been applied for a wide range of material behavior. The developed models can be classified as related to the single internal variable formalism based on the thermodynamic approach by Coleman and Gurtin. The extension of this theory to so-called dual internal variables opens up new avenues for the constitutive modeling of macroscopic material behavior. This paper reveals the distinction between constitutive modeling with single and dual internal variables using examples of heat conduction in rigid solids, linear thermoelasticity, and viscous fluids. A thermodynamically consistent framework for treating internal variables with as little a priori knowledge as possible is presented. This framework is based on the exploitation of the Clausius–Duhem inequality. Since the considered internal variables are “observable but not controllable”, only the Onsagerian procedure with the use of the extra entropy flux is appropriate for the derivation of evolution equations for internal variables. The key distinctions between single and dual internal variables are that the evolution equations are parabolic in the case of a single internal variable and hyperbolic if dual internal variables are employed.


Introduction
Constitutive modeling provides material-specific relationships between field quantities complementing balance laws. The emphasis is on developing mathematical models that can describe specific physical phenomena. The effectiveness of constitutive models is critical to the predictive capabilities of material computational models. Classical constitutive models, such as the Fourier law for heat flux and the Hooke law relating stress and strain, provide the description of homogeneous materials. Real materials have a microstructure composed of inhomogeneities of varying sizes, distributions, and properties. Even if the length scale of the microstructure is much smaller than the length scale of a body, the effect of the microstructure on overall body response may not be negligible. As a result, the influence of the microstructure should be considered in the macroscopic material description. Various approaches are used to describe the macroscopic behavior influenced by the internal structure of the materials. Homogenization methods result in the construction of effective media, which take into account the effects of the internal structure, e.g., [1][2][3][4][5][6]. More sophisticated generalized continua theories necessitate a large number of additional material parameters [7][8][9][10][11].
For a wide range of material behavior, phenomenological constitutive models with internal variables have been applied [12][13][14][15][16][17][18][19]. Even if multiple internal variables are introduced, the developed models can be classified as related to the single internal variable formalism based on the thermodynamic approach by Coleman and Gurtin [20]. The extension of this theory to so-called dual internal variables [21] opens up new avenues for the constitutive modeling of macroscopic material behavior.
The paper reveals the distinction between constitutive modeling with single and dual internal variables using examples of heat conduction in rigid solids, linear thermoelasticity, and viscous fluids. The primary challenge in implementing internal variable theory is the formulation of evolution equations for internal variables. The derivation of evolution equations is based on the exploitation of the dissipation inequality. Two common methods for exploiting the dissipation inequality are the Coleman-Noll procedure [22] and the Liu procedure [23]. Both procedures assume that the time rates of the state variables can be chosen arbitrarily. This assumption is not applicable for internal variables of state because they are "observable but not controllable" [13,24], which means that in contrast to dynamic degrees of freedom, internal variables cannot be governed by boundary conditions. Only the Onsagerian procedure with the use of thermodynamic forces and flux relationships is appropriate in this case [25].
We begin with the formulation of local balance laws and move on to the classic constitutive relations. Then, for heat conduction in rigid solids, linear thermoelasticity, and viscous fluids, a single internal variable framework is represented sequentially. The same phenomena are considered in the case of dual internal variables.
The main result of the paper is the demonstration of the ability of the dual internal variable approach to generate a hyperbolic evolution equation for internal variables that is applicable to various physical situations. It is shown in the case of heat conduction that a wave-like heat propagation can be described in a thermodynamically consistent manner. The causality issue for second gradient elasticity is avoided in the case of thermoelasticity. The Maxwell fluids model is recovered in the single internal variable case for viscous fluids, and extended equations of fluid motion are obtained using the dual internal variable approach.

Balance Laws
The behavior of a continuous medium is governed by balance laws. The local balance laws can be represented in the Eulerian (spatial) representation as follows (cf. [26]): • Conservation of mass: where ρ is the density field, v is the velocity field, and subscript t denotes the material time derivative; • Balance of linear momentum: where σ is the Cauchy stress tensor and f is a body force; • Balance of angular momentum: where upper index T denotes transposition; • Balance of energy: where e is the internal energy density, D is the strain rate tensor, q is heat flux, and h is a heat body source.
The entropy inequality supplements the balance laws where η is the entropy density field and θ is the absolute temperature.
In what follows, we consider the simplest possible situations: • Small deviation of temperature from the reference value; • Small-strain approximation; • Zero body sources; • Quadratic free energy.

Classical Constitutive Relations
As a starting point, we recall the classical constitutive models.

Heat Conduction in Rigid Solids
Because the matter density in rigid solids is constant by definition, dealing with internal energy per unit volume E = ρe and entropy per unit volume S = ρη is more convenient. The free energy per unit volume W = E − θS in rigid solids is only temperaturedependent, e.g., [27]: In the absence of body sources, the local balance law for energy in terms of free energy reads The second law of thermodynamics is written as Multiplying the second law (8) by temperature: and comparing the balance of energy (7) with Equation (9), we arrive at the inequality Here, "·" denotes the scalar product of vectors. Due to the arbitrariness of the temperature time rate, entropy per unit volume is defined as and the dissipation inequality (10) is reduced to This inequality is satisfied by the Fourier law: with the thermal conductivity k considered here as a positive constant. If the free energy per unit volume is a quadratic function where c p is the heat capacity, then the heat conduction equation for small deviations from the reference temperature θ 0 takes the classical form Classical heat conduction Equation (15) is valid for homogeneous isotropic materials obeying the Fourier law.

Linear Thermoelasticity
In the case of linear thermoelasticity, the free energy per unit volume depends on the strain tensor ε, the temperature, and, possibly, on the temperature gradient: Rewriting the balance of energy in terms of free energy: and multiplying the second law by temperature once more: we obtain the Clausius-Duhem inequality (no body sources) [28]: Here, ε t = D due to kinematic compatibility conditions. The substitution of the time derivative of the free energy into dissipation inequality (19) results in The arbitrariness of the time rates for strain, temperature, and its gradient yields This is the procedure used by Coleman and Noll [22] to exploit the Clausius-Duhem inequality. It leads to the definitions of entropy and stress in terms of free energy, and to the free energy independence of the temperature gradient. The rest of the dissipation inequality again, takes Fourier's law (13). For isotropic linear thermoelasticity, the free energy per unit volume is the quadratic function with Lamé coefficients λ and µ and the stress-temperature modulus β. For the stress, we have the relationship The equation of motion in terms of displacement u is then [28] u tt = µ∆u + (λ + µ)∇(divu) + β∇θ.
The balance of energy results in the heat conduction equation Thus, the two governing equations for linear thermoelasticity are coupled.

Linear Viscous Fluids
In the case of viscous fluids, we consider the isothermal situation. Here, the free energy density ψ = e − ηθ depends on the matter density ρ and the strain rate tensor D: The Clausius-Duhem inequality in the absence of body sources has the form [28] −ρψ t + σ : Applying the time derivative of the free energy density the decomposition of the stress tensor and mass conservation (1), we transform the Clausius-Duhem inequality into −ρ ∂ψ ∂D Here, p is pressure and σ v is the viscous stress tensor. The satisfaction of the dissipation inequality leads to the definition of pressure and to the free energy density independence of the strain rate tensor The rest of the Clausius-Duhem inequality determines viscous dissipation: For a compressible, linearly viscous (Newtonian) fluid the stress tensor is represented as [28] with the dilatational viscosity λ * and the shear viscosity µ * . Equations of motion follow from balances of mass and momentum (in the absence of body sources) [29]: These equations are generally referred to as the compressible Navier-Stokes equations.
The following steps are common in the development of classical constitutive models: • Choice of thermodynamic state space; • Calculation of the time derivative of free energy; • Combination of the first and second laws of thermodynamics to obtain the dissipation inequality; • Linear solution of the dissipation inequality.
This procedure will be expanded to include the influence of internal variables.

Simple Examples of Extended Constitutive Models
The three examples above all have one feature in common: they do not "feel" the internal structure of a material. For the description of more complex situations, constitutive models have to be extended. We indicate the simplest examples of such extensions.

Heat Conduction in Rigid Solids
The main goal of the Fourier law extension is to eliminate infinite speed for thermal disturbance propagation. Among the numerous non-Fourier models of heat conduction [30,31], the most well-known examples are the Cattaneo-Vernotte equation: and the Guyer-Krumhansl equation: where τ denotes a relaxation time and a is a material parameter. Despite the use of numerous alternative models and extensive research, the question of a suitable non-Fourier heat conduction model still remains open [32,33].

Linear Thermoelasticity
Generalized continuum theories, such as Mindlin's microelasticity [7] and Eringen's polar continua [34], extend elastic constitutive models completely and can be combined with thermal models [35]. However, even for the highly reduced so-called relaxed micromorphic theory [10], the large number of additional material parameters prevents their practical implementation. This is why strain gradient theories have gained popularity [36,37]. The majority of them, however, run into causality issues [38].
The Aifantis strain gradient model [39], in particular, is distinguished by the use of the strain's Laplacian to extend the stress-strain relationship: where c is an appropriate coefficient. This model is widely accepted in various applications [37,40,41].

Linear Viscous Fluids
The primary issue in hydrodynamics is turbulence [42]. There is currently no decent macroscopic description of turbulent motion [43]. Nonetheless, there are numerous rheological models that have no relationship to turbulence. A Maxwell fluid is a typical example of such a model. In this model, the stress tensor (36) is modified in the same way as in the Cattaneo-Vernotte equation [44,45] τσ t + σ = −pI + 2µ * D + λ * (trD)I. (42) Regardless of their differences, the constitutive models outlined above can all be considered as results of an internal variable approach. This is shown in the subsequent sections.

Single Internal Variable Framework
We consider an internal variable that is assumed to be capable of accounting for the influence of the internal structure to extend the constitutive models. The nature of the internal variable is not specified in advance. We investigate the impact of including the internal variable in constitutive models.
The procedure outlined in Section 3 also involves the following steps: • The derivation of evolution equations for internal variables as the solution of the dissipation inequality; • The possible elimination of an internal variable; • The analysis of the consequences for quadratic free energy and the linearization of results in particular cases.

Heat Conduction in Rigid Solids
We begin with a weakly nonlocal internal variable theory that comprises an internal variable α and its gradient into the thermodynamic state cf. [46,47]. This means that the free energy per unit volume W is specified as a function of the absolute temperature, the vectorial internal variable α, and its gradient [48]. With a quadratic free energy, we have where B and C are material parameters. In terms of free energy, the energy balance can be represented as where body forces are neglected for simplicity. Due to the dependence of the free energy on the internal variable, its time derivative is calculated using the chain rule As a result, the energy balance is rewritten as follows:

Dissipation Inequality
Keeping the presence of the internal variable in mind, we can write the second law of thermodynamics as where the entropy flux contains an extra entropy flux K cf. [12,[48][49][50][51]. In most cases, the extra entropy flux disappears, but this is not a necessary condition. We can rewrite the dissipation inequality by multiplying the second law (47) by temperature and taking into account energy balance Equation (46): If the extra entropy flux is chosen following the concept proposed by Maugin [12]: then the divergence term in the dissipation inequality disappears: The left-hand side of the obtained dissipation inequality now consists of the products of thermodynamic fluxes and forces.

Evolution Equation for the Single Internal Variable
The standard thermodynamic approach [52] provides the solution to the dissipation inequality by representing thermodynamic fluxes α t and (q + C∇α : α t ) as linear functions of conjugated thermodynamic forces: with and components M ij of the matrix M are considered as constants for simplicity. To ensure the non-negativity of entropy production, the symmetric part of the matrix M must be positive semidefinite [53,54]: As a result, the evolution equation for the internal variable α has the following form: Accordingly, we have the generalized heat flux of the relationship: The heat conduction equation is obtained in the framework of the single internal variable approach by replacing the heat flux in energy conservation Equation (46) with Equation (55): The right-hand side of the heat conduction equation depends on the internal variable. It is clear that the evolution equation of the internal variable (54) and heat conduction Equation (56) are coupled nonlinear parabolic equations.

Heat Flux as Internal Variable
It is worth noting that by using heat flux as the internal variable, we obtain the following for its evolution from Equation (54): It can be represented as an Guyer-Krumhansl-type equation: When C = 0, which means that the free energy is independent of the internal variable's gradient, the evolution equation for heat flux is reduced to a Cattaneo-Vernotte-type equation: The Green-Naghdi-type equation for heat flux is reached if, on the other hand, the free energy depends only on the gradient of the internal variable and is independent of the internal variable itself (i.e., B = 0): Thus, the Cattaneo-Vernotte equation and Guyer-Krumhansl equation are revealed in the framework of the single internal variable approach. The use of heat flux as an independent variable was studied in detail by Coleman et al. [55,56]. It was shown there that the Cattaneo telegrapher equation can be derived if, and only if, the internal energy does not depend on heat flux.

Linear Thermoelasticity
Expanding the thermoelastic state space using the tensorial internal variable α and its gradient, we can represent the free energy per unit volume W as a general sufficiently regular function: The balance of energy in terms of free energy keeps its form: and the second law multiplied by temperature is modified by the extra entropy flux The combination of Equations (62) and (63) results in the Clausius-Duhem inequality in the absence of body sources: As previously stated, the time derivative of the free energy should be substituted into Equation (64), yielding Again, the arbitrariness of the time rates for strain and temperature defines the stress and entropy: The rest of the dissipation inequality can be rearranged to The extra entropy flux is selected to eliminate the divergence term, following [12]: Then, the left-hand side of the dissipation inequality is transformed into the sum of the products of thermodynamic fluxes and forces:

Evolution Equation for Internal Variable
As before, the solution of the dissipation inequality is achieved by representing thermodynamic fluxes α t and (q + θK) as linear functions of conjugated thermodynamic forces. However, due to the difference in the tensorial order, both terms of Equation (71) should be non-negative separately. Therefore, the evolution equation for the internal variable α is For the heat flux, we have, correspondingly, Coefficients M i are considered as constants for simplicity.

Quadratic Free Energy: Isothermal Case
In the isothermal case, the free energy per unit volume is expressed as a simplified quadratic expression due to the complexity of including the third order tensor ∇α in full: Accordingly, the derivatives of the free energy with respect to the internal variable and its gradient are In this case, the evolution equation for the internal variable is a reaction-diffusion-like equation: which can be found under different names in numerous applications. This type of evolution equation is commonly used in phase field theory, e.g., [57]. In the nondissipative case, we have [58] C∆α − Aε − Bα = 0.
By definition, the stress tensor is the derivative of the free energy with respect to strain which determines the equation of motion as The equation of motion in terms of the displacement u is then It is possible to eliminate the internal variable from the equation of motion. First, we express divα in terms of the displacement Accounting for Equation (77), we obtain It follows from Equation (81) that Therefore, relationship (82) is rearranged to Finally, combining the latter equations with (81), we get It is this equation that has been analyzed in studies of dispersive wave propagation in microstructured solids [59][60][61][62]. It should be noted that dispersive wave Equation (85) confronts the causality problem [38].

Viscous Fluids
The Navier-Stokes equations describe the motion of homogeneous isotropic fluids. We use the internal variable formalism to take into account more complicated cases. In the single internal variable theory, the free energy density depends on the density ρ, temperature θ, and the internal variable α: The internal variable is not predetermined in advance. However, using the dissipation inequality (29), it is possible to derive an evolution equation of the internal variable in a general form.
The material time derivative of the free energy density is calculated using the chain rule Clausius-Duhem inequality (29) is represented as follows: The material time derivative of the density together with pressure relation (33), transforms dissipation inequality (88) into the linear combination of products of thermodynamic fluxes and forces: The evolution equation for the internal variable is derived using the obtained dissipation inequality.

Isothermal Case
We focus on the isothermal example to keep things as simple as possible. Dissipation inequality (90) in the isothermal case is reduced to As before, the linear solution of dissipation inequality (91) is given by i.e., where components M ij of the matrix M obey the same conditions of positive semidefiniteness as before. Relationships (93) and (94) establish the evolution equation for the internal variable and the thermodynamically consistent constitutive relation, respectively. These relations are not predetermined a priori or phenomenologically, but are rather the consequence of the dissipation inequality in the framework of the single internal variable theory.
In the case of a quadratic dependence of the free energy density the contribution of the internal variable on the right-hand side of Equations (93) and (94) is

Elimination of Internal Variable
There is no need to explicitly specify the internal variable because it can be eliminated from the dissipation inequality solution. After rearranging Equation (94) with (93), we have The material time derivative of relation (93) yields By substituting the time derivative of internal variable (97) into relationship (98), after some algebra, we obtain The influence of the internal variable is manifested in the structure of the relationship (99) and in the value of coefficient B in the quadratic dependency of the free energy density (95).
We notice that when we represent each tensor as the sum of hydrostatic and deviatoric components, their contributions to the dissipation inequality are orthogonal since A h : B d = 0 for arbitrary symmetric tensors. As a result, the solution of dissipation inequality (91) should be taken separately for the hydrostatic component (with the upper index h) and for the deviatoric component (marked by the upper index d) In the hydrostatic case, removing the internal variable from the solution of the dissipation inequality yields Similarly, the removal of the internal variable in the deviatoric example results in It is simple to observe that relationships (104) and (105) are simplified in the limit case and The combination of the last two equations: can be represented as constitutive Equation (36) for the Cauchy stress after identifying the values of coefficients since 3D h = (trD)I by definition. The Navier-Stokes equations clearly correspond to the considered limiting case.

Maxwell Fluids
It is now possible to observe how the Navier-Stokes equations can be expanded with the inclusion of the internal variable. To make matters easier, we set Aggregating Equations (104) and (105), we have We can represent this relationship as we arrive at a generalization of the constitutive equation for Maxwell fluids. The comparison with early models [44,45] and with recent research [63][64][65][66][67] shows that Equation (112) has a more complex structure. Any invariant material time derivative can be employed in Equation (112) [68].

Equations of Motion
The computation of the divergence of the stress tensor is required for the linear momentum balance ρv t = div σ.
Inserting the expressions of divergences into Equation (114), we obtain the equation which can be compared with a hyperbolic perturbation of the Navier-Stokes equations (c.f. [69]). As one can see, nonlinear terms appear due to the noncommutativity of the material time derivative and the divergence operators.
The classical constitutive relation holds in zero approximation In the first approximation, we have The corresponding divergences are as follows:

Incompressible Fluid
Consider the case of incompressible fluids as an example. The requirement of incom- reduces the relationships for divergences to Applying the identity we can represent the divergence of the first approximation of the stress tensor as Appeared linearized equation of motion can be compared to the gradient model by Aifantis [37] and the hyperstress model [70] ρv t = −∇p + µ * ∆v − ζ∆∆v.
Here, l is an internal length scale, and ζ is the difference between two constitutive moduli introduced in addition to classical viscosity µ * . Once more, a more complex equation of motion is produced by using internal variables. Because they result from the solution of the dissipation inequality, the obtained extended models of fluid motion are thermodynamically consistent. Generally, the internal variable is introduced to take into account the influence of motions at a microscale that cannot be accounted for explicitly.

Dual Internal Variables
As is demonstrated, the single internal variable framework validates the thermodynamic consistency of known constitutive models for heat conduction, linear thermoelasticity, and viscous fluids. Internal variables can sometimes be removed from equations of motion. The obtained constitutive models are constrained by the parabolicity of evolution equations of internal variables. To take things a step further, the dual internal variable theory should be applied.

Heat Conduction in Rigid Solids
The dual internal variable theory generalizes the single internal variable theory [21,62,71]. This is specified in [72,73] for heat conduction in solids. The free energy per unit volume is assumed to be dependent on temperature, internal variables α, β, and their space gradients in this approach.

Quadratic Free Energy
A quadratic function of its arguments is used to specify the free energy per unit volume: Then, the time derivative of the free energy is calculated as The obtained expression for the time rate of the free energy is then used to rewrite the equation for the energy balance (Equation (44)) The latter expression can be rearranged to The dissipation inequality in this case takes the form The extra entropy flux is chosen in order to eliminate the divergence term in the dissipation inequality [12]: Due to this choice, the dissipation inequality reduces to the sum of products of thermodynamic fluxes and forces:

Evolution Equations for Internal Variables
A linear solution of the dissipation inequality is constructed similarly to the case of single internal variable: where matrix L is assumed to be positive semidefinite as mentioned above. The evolution equations for internal variables α and β are represented as and the modified heat flux The heat conduction equation reads Nonlinear Equations (139)-(142) make up the closed system of equations to calculate the evolution of internal variables and temperature. In linear approximation, one of the internal variables (for example, β) can be eliminated.

Elimination of One Internal Variable
To obtain a single equation for the internal variable α, we compute the time derivative of evolution Equation (139) ignoring the nonlinear contributions for small temperature deviations from its reference value θ 0 : Next, the time rate of β can be expressed in terms of α by means of Equations (139) and (140): and the time rate of the Laplacian β: where the following notation is used: By substituting the relationships for the time derivatives of β and its Laplacian in the expression for the second time derivative of α, the single evolution equation for the internal variable α is obtained: The derived evolution equation is needlessly complicated. We consider its simplified version with vanishing values of coefficients B and F in the free energy dependence to reduce the complexity of the obtained evolution equation: Furthermore, if internal dissipation is negligible, i.e., L 13 = 0, L 11 = L 22 = 0, we arrive at It is a hyperbolic equation for the evolution of internal variable α if (− L 21 CD) is positive.
In the linear approximation, the heat conduction equation can be rewritten as Equations (148) and (149) are coupled. Together, they describe the heat conduction process in microstructured bodies, for example. In this case, the internal variable α represents the microtemperature, i.e., the temperature fluctuation around the macrotemperature. While the heat conduction equation for the macrotemperature remains parabolic, the evolution equation for the internal variable α is hyperbolic in contrast to the case of the single internal variable. This is the key distinction between frameworks with single and dual internal variables.

Linear Thermoelasticity
For linear thermoelasticity, the internal variable theory is easily generalized to the situation where there are two internal variables. Consider the free energy per unit volume W as a function of two internal variables, α and β, each of which is a second-order tensor: θ, α, ∇α, β, ∇β).
The energy balance in terms of free energy preserves its structure: and the second law multiplied by temperature is modified by the extra entropy flux: The combination of Equations (151) and (152) results in the Clausius-Duhem inequality in the absence of body sources As previously, the time derivative of the free energy should be substituted into Equation (153), yielding Again, the arbitrariness of the time rates for strain and temperature defines the stress and entropy: The rest of the dissipation inequality (157) can be rearranged to The extra entropy flux is selected to eliminate the divergence term, following [12]: Then, the left-hand side of the dissipation inequality is transformed to the sum of products of thermodynamic fluxes and forces:

Dissipation Inequality: Isothermal Case
Until now, the introduction of a dual internal variable looks like a repetition of the single internal variable formalism. However, the dissipation inequality is more general even in the isothermal case and its linear solution where components M 11 , . . . , M 22 of the linear operator M are dependent on state variables [74]. As previously stated, the matrix M should be positive semidefinite.

Quadratic Free Energy
Representing the free energy per unit volume as the quadratic function we obtain four derivatives of the free energy The evolution equations for internal variables are specified as As before, it is possible to derive a single evolution equation for the internal variable α: A simplified version of this evolution equation is obtained for B = 0 and F = 0: In the nondissipative case (M 11 = M 22 = 0 and M 12 = −M 21 ), the evolution equation is reduced to a hyperbolic equation: By definition, the stress tensor is the derivative of the free energy with respect to strain: which determined the equation of motion as The equation of motion in terms of displacement u is then This is a hyperbolic equation of higher order as is discussed in [58][59][60][61][62]71,[75][76][77]. It is obvious that Equation (172) has no issues with causality and is valid for the propagation of dispersive waves [77]. In contrast, the causality issue is unavoidable in the single internal variable approach.

Viscous Fluids
Now, we study the difference between single and dual internal variable descriptions in the case of viscous fluids. In the dual internal variable approach, it is supposed that the free energy density depends on internal variables α and β: θ, α, β).
We compute the time derivative of the free energy: to use it in the Clausius-Duhem inequality: If the free energy density is quadratic ψ(. . . , α, ∇β) = . . .
then the dissipation inequality is transformed into the sum of products of thermodynamic fluxes and forces taking into account Equation (33): As before, the representation of thermodynamic fluxes as linear functions of conjugated thermodynamic forces leads to the solution of the dissipation inequality (177).

Isothermal Case
The dissipation inequality is even simpler in the isothermal situation The linear solution to the dissipation inequality is obtained in the same manner as in the case of a single internal variable: where with the same condition concerning the positive semidefiniteness of the matrix L as for the single internal variable case. As a result, the constitutive relationship for stress is depicted as Evolution equations for internal variables α and β take the form Both the constitutive relationship (181) and the evolution equations are coupled. The internal variables α and β are dual, as shown by Equations (182) and (183).

Elimination of One Internal Variable
We can eliminate one of the two internal variables. We take the time differentiation of Equation (182) to obtain and then develop expressions for β t in terms of α using the relation which follows from Equation (182). Thus, we have We obtain the evolution equation of the internal variable α in its own terms by substituting equality (186) into Equation (184): The internal variable α alone can also be used to rewrite constitutive relation (181): The last equation can be resolved for α t : which gives another expression for the evolution equation of the internal variable α: When we compare Equations (187) and (190), we get the solution to the dissipation inequality in the form In the dual internal variable approach, Equations (187) and (191) provide the evolution equation for the internal variable α and the constitutive relation for fluid flow.

Asymptotics
It is difficult to analyze the obtained constitutive relationship in its entirety. However, the relationship is divided into two parts that can be considered separately. In fact, the coefficient D describes the impact of the second internal variable on the overall behavior. If this coefficient is small, than the right-hand side of Equation (191) can be neglected, which results in The evolution equation for the internal variable α is then reduced to and the constitutive relation of Equation (188) is represented as We can eliminate the internal variable α by combining the obtained relationships which produces a variant of the constitutive equation for the case of a single internal variable, as expected.
In the opposite case, when the coefficient D is large, we obtain a different constitutive equation: The equation of motion for large values of D is This equation resembles the Navier-Stokes equation of motion. The last term on the righthand side, however, is dependent on the divergence of the internal variable. The situation reduces to the classical incompressible case in the absence of internal variables. The evolution equation for the internal variable α (187) is based on Equation (185), which can be hyperbolic if the internal variable β is interpreted as the gradient ∇α: and coefficients L 23 and D have opposite signs.

Conclusions
Constitutive modeling with internal variables is based on the premise that traditional constitutive models are merely a first approximation in describing actual phenomena and should be enhanced for broader applicability. Internal variables represent what we are unable to directly describe phenomenologically. The widely exploited phase-field approaches [78][79][80] are, in fact, examples of the use of internal variables [57,81]. The same is true for the gradient elasticity theory [58]. Moreover, the Mindlin microelasticity and Cosserat microrotation theories can be interpreted in terms of dual internal variables [71]. The duality of internal variables refers to the use of two equally valuable coupled internal variables to describe the same phenomenon rather than the use of multiple internal variables to describe different phenomena [21].
The paper presents a thermodynamically consistent framework for treating internal variables with as little a priori knowledge as possible. This framework is based on the exploitation of the Clausius-Duhem inequality. Because the considered internal variables are "observable but not controllable", only the Onsagerian procedure with the use of the extra entropy flux is appropriate for deriving of evolution equations for internal variables.
The most noticeable effect is the clear distinction between constitutive models with single and dual internal variables. It is shown that the evolution equations for internal variables are parabolic in the case of a single internal variable and may be hyperbolic if dual internal variables are employed. This conclusion holds true for all considered examples.
In the case of heat conduction in rigid solids, the constitutive model with a single internal variable provides both the Cattaneo-Vernotte and the Guyer-Krumhansl equations if the internal variable is interpreted as heat flux. However, both the heat conduction equation and the internal variable evolution equation are parabolic. Only in the case of dual internal variables is it possible to obtain a hyperbolic evolution equation for the internal variable after one of them is eliminated. The remaining internal variable can be interpreted as a microtemperature (the fluctuation around mean temperature). Although the heat conduction equation remains parabolic, the coupling with the microtemperature evolution causes wave-like heat propagation [75,82,83].
The constitutive models for thermoelasticity are illustrated using dispersive wave propagation as an example. The internal variables here can be interpreted as microdeformations [62,71,75]. In the case of a single internal variable, the corresponding models confront the issue of causality [38]. This issue can only be eliminated by employing dual internal variables [77].
The extension of the Navier-Stokes equations yields constitutive equations for Maxwell fluids in the framework of a single internal variable [84]. In contrast, the dual internal variable approach is able to produce Reynolds-type equations of motion [43] with a possibly hyperbolic evolution equation for an internal variable.
It should be emphasized that the constitutive models with internal variables are typically nonlinear and generalize existing models [85][86][87][88]. This opens up a wide range of possibilities for applying internal variable theory to constitutive modeling of complex phenomena. The paper includes simplified examples of such models so that they can be compared to known models.